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ABSTRACT 

We develop an analytical model for the accretion and gravitational drag on a point mass that moves hyper- 
sonically in the midplane of a gaseous disk with a Gaussian vertical density stratification. Such a model is 
of interest for studying the interaction between a planet and a protoplanetary disk, as well as the dynamical 
decay of massive black holes in galactic nuclei. The model considers that the flow is ballistic, and gives fully 
analytical expressions for both the accretion rate onto the point mass, and the gravitational drag it suffers. The 
expressions are further simplified by taking the limits of a thick, and of a thin disk. The results for the thick 
disk reduce correctly to those for a uniform density environment (Canto et al. 2011). We find that for a thin 
disk (small vertical scaleheight compared to the gravitational radius) the accretion rate is proportional to the 
mass of the moving object and to the surface density of the disk, while the drag force is independent of the 
velocity of the object. The gravitational deceleration of the hypersonic perturber in a thin disk was found to be 
independent of its parameters (i.e. mass or velocity) and depends only on the surface mass density of the disk. 
The predictions of the model are compared to the results of three-dimensional hydrodynamical simulations, 
with a reasonable agreement. 

Subject headings: black hole physics — hydrodynamics — stars: formation — ISM: clouds — ISM: kinematics 
and dynamics 



1. INTRODUCTION 

A gravitational point-like particle moving in a background 
gas can capture ambient matter. In addition, the particle may 
experience a drag due to the interaction with the wake that it 
induces in the medium. Both accretion and dynamical fric- 
tion are ubiquitous phenomena in astrophysics, from proto- 
stars embedded in molecular clumps to supermassive black 
holes at the center of galaxies (e.g., Sanchez-Salcedo 2012). 

Bondi, Hoyle and Lyttleton (Hoyle & Lyttleton 1939, 
1940a, 1940b, 1940c; Bondi 1952) treated the problem by 
assuming a point mass moving at constant speed within a uni- 
form gaseous media. By neglecting the pressure of gas ex- 
cept on the downstream axis, Bondi & Hoyle (1944) derived 
a formula that gives the accretion rate as a function of the 
Mach number for highly supersonic flows. Later on, Bondi 
(1952) suggested an approximate formula for any Mach num- 
ber. Some authors have used Bondi's formula to account for 
the initial mass function of stars as a competition for gas ac- 
cretion by protostars (Bonnell et al. 2001 a,b; Klessen & Burk- 
ert 2000, 2001) and to estimate the rate at which supermassive 
black holes accrete mass (e.g., Di Matteo et al. 2001). 

The loss of momentum of the object relative to the cloud 
due to its gravitational interaction with its wake was also en- 
visaged in Bondi & Hoyle (1944) using their line-accretion 
model. Traditionally, however, the gravitational drag force 
has been derived by calculating the density structure of the 
wake in linear perturbation theory (Dokuchaev 1964; Ruder- 
man & Spiegel 1971; Just & Kegel 1990; Ostriker 1999; Kim 
& Kim 2007; Sanchez-Salcedo 2009, 2012; Namouni 2010). 
In this approach, a minimum radius r m i n must be introduced 
because the linear approximation is not valid close enough to 
the perturber. Recently, Canto et al. (2011) were able to cal- 
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culate the contribution of the nonlinear inner part of the wake 
to the gravitational drag on hypersonic perturbers by using the 
ballistic orbit theory. 

All the abovementioned papers have considered that the 
surrounding gaseous medium is initially uniform, of constant 
density pq. In this case and for a poly tropic gas, the accre- 
tion rate and the gravitational drag are largely characterized 
by the Mach number. Because of the assumption that the 
medium is infinite, the gravitational drag on supersonic bodies 
increases logarithmically in time (e.g., Ostriker 1999). In real 
life, however, systems are finite in size. Whilst the problem 
of the dynamical friction in non-homogeneous (usually flat- 
tened) collisionless systems has received considerable atten- 
tion (eg., Binney 1977; Tremaine & Weinberg 1984; Maoz E. 
1993; Colpi et al. 1999; Just & Penarrubia 2005), less studied 
is the case for gaseous systems. Sanchez-Salcedo & Branden- 
brug (2001) simulated the orbital decay of a satellite in orbit 
around a Plummer sphere of gas and found that the "local ap- 
proximation", that is estimating the drag force at the present 
location of the perturber as if the medium were homogeneous 
(but taking appropriately the Coulomb logarithm) is very suc- 
cessful. 

In this paper we study the gravitational drag on a body 
moving in the midplane of a vertically-stratified gaseous disk. 
It has been recently noted that the dynamical friction in a 
gaseous slab could be useful to describe the interaction be- 
tween a planet in an eccentric orbit and the protoplanetary 
disk (Muto et al. 2011). These authors argue that the stan- 
dard analysis of the interaction between a planet and the disk, 
which assumes that the planet is corotating with the disk, is 
not valid when the planet eccentricity exceeds the disk as- 
pect ratio. Instead, for moderate and large eccentricities, they 
suggest it is more convenient to use the dynamical friction 
formula. They derived the dynamical friction force felt by a 
perturber moving in an infinitesimally thin disk of constant 
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surface density and calculated the migration and eccentricity 
damping timescales. When the dynamical friction formula is 
applied to planets with moderate eccentricities, the migration 
and eccentricity damping timescales are in agreement with 
those found in previous works. 

Gaseous dynamical friction is also relevant to study the 
disk-planet interaction in highly inclined systems. Rein 
(2012) computed the time-scales associated with migration 
and inclination damping of planets due to dynamical friction 
with the protoplanetary disk in the linear approximation and 
compared the results with simulations. In order to model the 
gravitational interaction, he uses a large softening radius, so 
that their simulations also lie in the linear regime as well. This 
is justified as long as the physical radius of the planet r p is 
much larger than its accretion radius (Ra = 2GM/v 2 , where 
v is the velocity of planet relative to the gas in the disk). For 
a hot Jupiter with a mass of 100 Earth masses and a size of 10 
Earth radii moving on an inclined circular orbit around a star, 
the condition r. p > Ra implies inclinations between the orbit 
of the planet and the plane of the disk larger than ~ 71°. In 
the present work, we are interested in modelling the nonlinear 
contribution of the gravitational wake to the drag force on an 
accreting point-mass moving in non-inclined orbits. 

A second area of application is in the orbital decay of mas- 
sive black holes residing in galactic nuclei. The formation 
of massive black hole binaries following galaxy mergers is a 
natural consequence of the hierarchical scenario (e.g., Begel- 
man et al. 1980). The gravitational torque due to dynamical 
friction with an accretion disk formed by gas that was fun- 
nelled to the central 100 pc, is expected to be responsible for 
the in-spiral toward the center of the merger remnant down to 
pc separations (Escala et al. 2005; Dotti et al. 2007; Cuadra 
et al. 2009; Khan et al. 2011; Preto et al. 2011; Roedig et 
al. 2012). Recently, Nixon et al. (2011a,b) showed that the 
binary orbital angular momentum and energy are more effi- 
ciently removed in a retrograde accretion disk than they are 
in a prograde accretion disk because there are no orbital res- 
onances between the binary and the disk. Since there is no 
compelling reason to assume that the binary and the disk ro- 
tation are initially parallel, it is quite possible that retrograde 
disks strongly promote coalescence of the black hole binary. 

So far, most of the semi-analytical models aimed to study 
the assembly of black holes use a prescription for the dynam- 
ical drag (e.g., Tanaka & Haiman 2009). Given the impor- 
tance of having analytical formulae to be applied for plan- 
ets or black holes in eccentric or counter-rotating orbits, the 
strategy in our paper is to provide analytical estimates for the 
gravitational drag on an accreting object embedded in a three- 
dimensional disk of gas and to check its validity through nu- 
merical simulations. We will focus on highly supersonic per- 
turbers. For instance, if a perturber is in orbit within a disk 
whose aspect ratio is constant with radius, the relative veloc- 
ity between the perturber and its surroundings will be always 
supersonic if its eccentricity exceeds the disk aspect ratio (see, 
for instance, Muto et al. 2011). Moreover, it is clear that per- 
turbers in retrograde disks also move supersonically, even if 
the orbit is circular. We will concentrate on a Gaussian ver- 
tical profile, but the method is general and can be applied for 
any decaying profile. 

2. THE FREE STREAMING MODEL 

We consider the problem of a point mass M moving hyper- 
sonically at velocity vq in a stratified, gaseous medium. As a 
result of the gravitational pull produced by M the otherwise 
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Fig. 1. — Schematic diagram of the geometry and reference frame used 
in the model (the mass M is at rest). Far from the source a fluid parcel has 
a velocity vq (parallel to the a'-axis), and an impact parameter £, but as it 
approaches the point mass its trajectory is deflected as illustrated in panel 
(a). In panel (b) the angle ip is shown in a view perpendicular to the relative 
motion of M with respect to the environment. 

straight trajectory of a fluid parcel is bent, and a wake behind 
the perturber is formed (see the top panel of Figure [TJ. In the 
hypersonic limit, the solution for the trajectory of a fluid par- 
cel (streamline) can be obtained neglecting the pressure force 
and considering that the motion is ballistic. By hypersonic 
we mean that the dynamical pressure is much larger than the 
thermal pressure, i.e. Vq ^> c 2 , where c s is the sound speed, 
so that the Mach number M. satisfies M 2 3> 1. The case of a 
uniform environment was treated in |Canto et al.| ( [20lT) and a 
comparison of the analytic model and a numerical simulation 
was presented. 

In the present paper we consider that the point mass moves 
in the mid-plane of a plane-parallel structure, whose density 
decays in the vertical direction y. For illustration, we will 
assume a rather realistic Gaussian vertical density profile of 
the form 



p(y) = p e-y 2 l h2 = Po e-^ cos ^ 2 / h2 



(1) 



where p is the midplane density, h is the scaleheight, £ is the 
fluid parcel impact parameter, and the angle ip is measured 
from the y-axis (see Figure fT). We adopt the same reference 
frame as in |Cant6 et al.| ( |201 1[ ), in which the point mass M 
is at rest. In this reference frame, the motion can be modeled 
as a streaming environment with a constant velocity vo far 
upstream of the source. The undisturbed, streaming velocity 
v lies in the positive x direction. 

In a hypersonic regime, the trajectories of the parcels, sub- 
ject to the gravitational field of M, have initially positive en- 
ergy (per unit mass) E = v 2 /2 and therefore they are hy- 
perbolic. The velocity of each parcel of gas is the same as 
Canto et al. ( 201 1[ ), that is, the solution of the hyperbolic 



m 



trajectory after imposing the upstream vq boundary condi- 
tion and requiring angular momentum conservation along the 
streamlines. Our model assumes that when the fluid encoun- 
ters the x-axis, the other velocity components (v y and v z ) are 
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thermalized at a shock, and the corresponding energy is im- 
mediately radiated away. This condition results in a decrease 
of the kinetic energy from Eq = Vq/2 at x — > — oo to 



E u 



1 



(2) 



where £p is the gravi tational radius ( Bisnovatyi-Kogan et al. 
[19791 |Cant6 et al.|201 1| >, defined as 

<b-^. (3) 

v o 

From equation d2]i it is clear that all the material with an im- 
pact parameter f < 2£o after the shock has a negative energy 
and thus remains gravitationally bound (i.e., it will eventually 
be accreted onto the point mass). 

From equation Q we can also see that the material that 
is gravitationally unbound when entering the wake will flow 
away along the a>axis reaching infinity with a velocity 
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2. 1 . The mass accretion rate and the gravitational drag 

The accretion rate can be calculated by integrating over all 
the material intercepted by a cylinder with an impact parame- 
ter up to 2£ : 

M m .= I 2^ / 5(£)f £d£, (5) 



where p(£) is the -0-averaged density for a given impact pa- 
rameter £. This average is obtained as: 

HO = jf" I* P&4>) W = Poe-« 2/2 ' l2 I , (6) 



where I n (a;) is the modified Bessel function of the first kind. 
Substituting equation (|6|l into equation |5]) and performing the 
integral one obtains: 

h 2 I \ h? 



M acc = 4-7T po v £ e 



(7) 

The drag force felt by the point mass can be calculated con- 
sidering the net loss of momentum flux along the x-direction. 
Material coming into the system from x — > — oo carries a pvo 
momentum flux, while material leaving the system at x — s- oo 
(fluid parcels with £ > 2£o) takes a pvoo momentum flux, thus 
the drag can be obtained as 



(8) 



Fd = V + / 27T V p(£) (Vq - Woo) £ d£. 



26) 



The first term in the right-hand side is the drag caused by ac- 
cretion of material and the second term is a purely gravita- 
tional drag by the material in the wake that will not be ac- 
creted. In order to obtain a closed expression for it is 
convenient to introduce the new variables I = £/(2fo)> and 
A = £o/ h. With these variables and replacing Eqs. Q and ^ 
into equation dSl the drag force can be written as: 



= M acc V 



8irp vl £ 2 x 



I (2A 2 I 2 ) I - (I 2 - 1) 



1/2' 



dl .(9) 



The expressions for the mass accretion rate (equation |7J, 
and for the gravitational drag (equation[9]l are easy to compute 
numerically. Moreover, they can be further simplified in the 
limiting cases of a thin and thick disk, as shown below. 



2.2. Thick disk limit: A< 1 

If the scaleheight of the disk is much larger than the gravita- 
tional radius (that is, A = (£,o/h) <C 1), the Io Bessel function 
in equation |7) is close to unity, whereas the Ii is close to 
zero. After substituting £o by its definition in Equation ([3]), 
the accretion rate is approximately given by: 



ft I, 



acc, thick 



4tt (GMfpo 



(10) 



which is th e same exp ression obtained for the uniform density 
case ( |Cant6 et al.|2011| >. 

The mass accretion rate can be written in terms of the disk 
surface density given by 



E = 



p{y) dy 



-y 2 /h 2 



dy=y / Trp h, (11) 



as : 



ftl 



acc. thick 



Vtt (GM) 2 £ 



hvl 



(12) 



We may obtain an analytical expression for the drag force 
by using the following approximations 



-lT A 



Io (2r 2 ) 



e~ 2r I„(2r 2 ) 



; 2r 2 < 1, 

; 2t 2 »i, 



(13a) 
(13b) 



to split the integral in equation (|9]l in two parts. A first integral 
over I in the range 1 < I < lo (such that 2A 2 I 2 < 1), and the 
second x-integral from xq to infinity (such that 2A 2 I 2 > 1), 
where 

/0 = 2^A (U) 



is the junction of the asymptotic expressions in Eqs. (13 1 
More specifically, 



-2\ z r 



I (2A 2 Z 2 ) I -(I 2 - 1) 



1/2' 



dl 



(I 2 
l- 



,1/2' 



I 2 -I 



dl 

,1/2' 



dl. 



(15) 



Performing the integrals, and using the accretion rate of equa- 
tion ( 10 1, the total drag becomes 



-Fd.thick ~ 47r po v 2 £ 2 



2 Iq arctan 



1 



In lo + \ Hi - 1 



(16) 



Remind that in terms of physical parameters, we have that 
£o = (GM) 2 /d 4 and l = hvl/(2^GM). The associated 
Coulomb logarithm, defined as In A = F&/ ' (Air po(; 2 v 2 ), cor- 
responds to the term in square brackets in Eq. ( fT*5j ). We see 
that the ambiguity in the definition of the minimum scale of 
interaction that appears in linear theory (e.g., Ostriker 1999 
for the homogeneous medium, or Rein 2012 for the case of a 
disk), is removed in this approach. Moreover, no additional 
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cut-off scale needs to be introduced at large distances. Us- 
ing the hydrodynamical approach in the linear regime, Rein 
(2012) identified a Coulomb logarithm but it is not the same 
as in the present work. 

A yet simpler expression can be obtained taking the limit of 
l > 1 (i.e. A < 1; the thick disk limit) : 



F 



d, thick 



4irp {GMf 



In 



(17) 



The error made by using the above expression instead of the 
exact formula (Eq. [9|, is less than 2% for A < 0.5. In terms 
of the disk surface density, the drag force in the thick disk 
approximation can be expressed as 



thick 



40FE(GA/) 2 



In 



e 3 / 2 h 



(18) 



where the factor e 3 / 2 comes from inserting the term 3/2 that 
appears in the right-hand side of Equation (17i in the argu- 
ment of the logarithm. 

For a homogeneous medium, Canto et al. (2011) demon- 
strated that the drag force for hypersonic bodies can be written 
as 



Fa 



4-kpq {GMf 



In 



(19) 



is the largest impa ct p arameter and 



where r 

Comparing equations 
in the case of a stratified medium. 



18 



and 



19 1 we find that, 



^ne ambiguity in the maxi- 
mum impact parameter is removed. Adopting r m i n = ^/e^Q/2 
and using Eq. ( 18 1, we find r max — e 2 h/ ~ 2.1 h for a 
Gaussian disk. Hence, this implies that the gravitational drag 
on a supersonic body that is dropped suddenly at t — in a 
stratified background, saturates asymptotically with time to a 
constant value. This also holds true if the disk is infinitesi- 
mally thin (see the next subsection and Muto et al. 201 1). 

2.3. Thin disk limit: X^$> 1 

For a thin disk (h <C or A » 1), one can use the 
approximation of equation ( |13a| i in equation (|7]i to obtain the 
mass accretion rate: 



Ma 



;,thii 



Ay/n^ohpo Vq, 



(20) 



or in terms of the disk surface density (equation 111, after 
substitution of equation Q, as 

4GMS 

:,thin « 4 £ £ = — ~ ■ (21) 



t'o 



Notice that the accretion rate has a much weaker dependence 
on the velocity of the perturber vo than in the thick disk or 
uniform density cases (Equation[T2|), and that the dependence 
on the disk parameters is nicely folded into a single parameter 
(i.e., the surface density). 

The net gravitational drag (including the drag due to ac- 
cretion) on a body moving in a thin disk can be estimated 
combining equations (9 13a and 20 1, yielding : 



fih£ Q = 2ttGM£. (22) 



Thus, the drag in a thin disk does not depend on v (i.e., the 
relative velocity between the point mass and the surrounding 
environment), but solely on the surface density £ of the disk 
and the mass M of the perturber. Furthermore, the gravita- 
tional deceleration on the hypersonic perturber F^/M is in- 
dependent of any of the parameters of the perturber (i.e. its 



mass or velocity) and depends, only, on the surface mass den- 
sity of the disk. 

It is common to compare the drag force in gaseous and col- 
lisionless media. Following Binney & Tremaine (1987), we 
have derived the dynamical friction force in a collisionless 
medium but, instead of assuming that the background par- 
ticles are distributed homogeneously in a three-dimensional 
medium, we assume that they are distributed in the (x, z)- 
plane, i.e. two-dimensional geometry. We find that, if the 
perturber moves at a velocity much larger than the velocity 
dispersion of background particles, the dynamical friction is 
also given by Eq. p2| ). We must stress here that the struc- 
ture of the wake and the underlying physics in gaseous media 
(where we may have gas accretion onto the perturber) and in 
collisionless media (no accretion occurs) are very different. 

3. NUMERICAL SIMULATIONS 
3.1. The setup 

We have used the adaptive grid code YGUAZU-A (Raga 
et al.||2000| to perform a set of numerical simulations to be 
compared with the analytical model presented above. The 
version of the code employed solves the isothermal hydrody- 
namic equations on a three-dimensional Cartesian grid. The 
equations are solved with a second order implementation of 
the "flux vector splitting" algorithm of van Leer ( 1982}. 

We use the same reference system as that described in Sj2| 
in which the point mass is at rest at the center of the sys- 
tem. In units of £o> tne computational domain covers a re- 
gion of [(—7.5, 7.5), (0, 7.5), (—7.5, 7.5)], in x, y and z, re- 
spectively. The relevant scales that need to be solved are 
£o and h, the gravitational radius and the scaleheight of the 
disk, respectively. The domain is discretized in a five-level 
binary adaptive grid with a maximum resolution equivalent 
to 512 x 256 x 512 cells in a uniform mesh. Simulations at 
half the resolution give similar results for the accretion rate 
and the drag, but with a smoother structure of the wake, due 
to numerical diffusion. The gravitational attraction due to a 
point mass centered at (0, 0, 0) is added as an external source 
term. To avoid numerical artifacts a softening length of 10 
cells is used for the gravitational force. The accretion onto 
the point mass is emulated by artificially keeping (at every 
timestep) a low density inside a hemisphere of radius 0.15 £o 
(~ 5 cells) centered around the point mass. 

Different parameters have been used in the literature to 
measure the degree of non-linearity A. In the case of the 
Bondi-Hoyle accretion problem, we will use A = 2^/r s , 
which is the ratio between the accretion radius (2£ ) and the 
softening radius r s . In our simulations r s = 0.3£o and hence 
A = 2£o/?'s — 6. In astronomical systems, A may range 
between 1 and 100 for Earth-like and Jupiter-like planets, up 
to much higher values for black holes. However, a value of 
A ~ 6 is large enough to test the success of our analytical 
formalism in a situation in which non-linearity is important. 

Since the problem has mirror symmetry on the y = plane, 
we only simulate one half of the wake (for positive y), impos- 
ing a reflective boundary condition at y = 0. The subsequent 
analysis is done considering the mirror symmetry. The do- 
main is initially filled with a streaming environment with a 
density given by the Gaussian profile of equation ([T]), and a 
constant velocity vq in the x direction. An inflow condition 
(with the appropriate density profile and velocity) is imposed 
on the x = — 7.5£o boundary in order to replenish the stream- 
ing environment. Outflow conditions are imposed on the re- 
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TABLE 1 
Parameters of the simulations 3 . 



Model 


c s 


Mach number 


h [Ho] 


Regime 


A 


0.2 


5 


10 5 


thick disk 


B 


0.2 


5 


2.0 


thick disk 


C 


0.2 


5 


0.5 


thin disk 


D 


0.1 


10 


0.5 


thin disk 



a In all models po = 1, and vo = 1.0. 

maining boundaries. 

Because the flow is isothermal, the density gradient trans- 
lates into a vertical pressure gradient that tends to destroy the 
disk, which is particularly problematic for the thin disk mod- 
els (which have a larger pressure gradient). To avoid this we 
have added a vertical external force (in the —y direction) that 
corresponds to the gravity force needed to maintain hydro- 
static equilibrium. Additionally, to avoid numerical problems 
due to the large range of density values needed for the Gaus- 
sian profile, we have set a lower limit beyond a height of 
y = 3£o- Past this height the density is homogeneous and 
the external vertical gravity is turned off. This has no no- 
ticeable consequences on the results as most of the mass and 
momentum lie closer to the midplane. 

The parameters of the four models are summarized in Ta- 
ble [T] In all the models, we take units such that p a = 1 and 
vq = 1. The first two simulations (A and B) correspond to a 
thick disk. Given the large scaleheight used in model A, the 
ambient density is close to constant inside the domain, thus 
the model is basically a 3D version (albeit at a lower reso- 
lution due to computational constraints) of the simulation in 
|Cant6 et al.| l |201l) . The other two models (C and D) corre- 
spond to a thin disk with the same scaleheight but different 
sound speeds. Models C and D have Mach numbers of 5 and 
10, respectively. These two models were computed to test the 
result (of the analytic model) that the drag is independent of 
the Mach number of the flow. 

3.2. Results 

We allow the four models to run up to an evolutionary 
time of t = 55£o/vo, ensuring that a quasi-stationary state is 
achieved. Density maps for all the models after this integra- 
tion time are displayed in Figure [2] In this figure we present a 
zoom of a region close to the point mass. The snapshots show 
the highly variable flow structure, formed downstream of the 
perturber centered at (0, 0, 0). 

In order to compare the numerical simulations with the free- 
streaming model, we have taken the density and velocity out- 
puts and computed the mass flux (per unit of time) M and 
the net a; -momentum flux, over cylinders of varying radii (£). 
More specifically, 



M(0 



*b(0 



pv ■ dS, 



pv x v ■ dS, 



(23) 



(24) 



where S% is the surface of a cylinder with axis along the x- 
axis, radius £ and whose caps are placed at the boundaries of 
our domain, at x = — 7.5£o and at x = 7.5£o- In a steady 
state, M should not depend on £ and must coincide with the 
accretion of mass onto the perturber M acc - On the other hand, 



lOi, (171, 



Fd(0 is the contribution of the mass inside a cylinder of ra- 
dius £ to the drag force. It includes both the force due to 
accretion of momentum over the body surface and the gravi- 
tational force on the perturber by its own wake. The drag force 
on the perturber Fa corresponds to when all the mass in 
the wake is included, that is, in the limit £ — > oo. The results 
can be compared with the predictions of equations 
|2T} and (^. 

The results for M as a function of the cylindrical radius 
(i.e., the impact parameter) are presented in Figure [3] This 
figure gives the average M from t = 51 to 55, in units of 
£o/ w o> with error bars denoting the maximum variability at 
each radius. The highly variable mass flux, sign of a very tur- 
bulent wake, evidenced by the large error bars, is consistent 



with previous calculations (see for instance Canto et al. 201 1[ ). 
The solid line in the figures correspond to the asymptotic ex- 
pressions for the thick (A, B) and thin (C, D) disk models, 
Equations ( fl0| ) and ( |20| , respectively. We have added the full 
solution for M (Equation^, which is only appreciably differ- 
ent from the solid line in Figure [3jb). We can also see from 
Figure [3] that M becomes more or less independent of the im- 
pact parameter after £ = 2 £0, the impact parameter beyond 
which the material remains unbound (see eq. [3}. 

It is also noticeable that there is a systematic difference to- 
wards smaller accretion rates compared to the analytical pre- 
diction. The reason for this offset is the limited length of the 
computational box. Material entering the downstream wake 
will travel a distance (see |Canto et al.|201 l| l 



2£o 



1 



(25) 



along the positive x-axis before turning around and falling to- 
wards the point mass. With the outflow boundary condition 
placed at x = 7.5 £0, all the mass that leaves the computa- 
tional box is lost. From equation ( 25 1 we can see that only 



material with £ < 1.77 £0 is captured, while all the mass with 
an impact parameter £ > 1.77 £0 leaves the domain and does 
not contribute to the mass accretion. The effect of the limited 
domain size can be quantified by integrating the mass accre- 
tion (equation|5]l only up to an impact parameter 1.77 £0. The 
result is shown as a dashed line in the plots of Figure[3] agree- 
ing well with the results from the simulations. 

The flux Fd as a function of impact parameter (averaged 
over the same time window as the mass flux, see above) is 
presented in Figure |4] The gravitational drag in our simula- 
tions also shows a saturation for £ > 2 £0, for all the models. 

Comparing the two thin disk models in Figs.[4jc,d) we can 
see that indeed Fa is basically independent of the Mach num- 
ber (provided that the flow is hypersonic). 

The solid lines shown in Fig. [4] correspond to the analytical 
approximations given in equations ( fl7] > and ( |22| , for the thick 
and thin disk models, respectively. We have also included the 
full solution of Equation |9]) as a dotted line, which overlaps 
the solid lines except for model B. Since the net drag Fa has a 
term M acc uo (see Eq.[9|, we need to take into account that the 
accretion rate is slightly smaller than predicted by the model 
because of the limited extent of the computational box in the 
a; -direction. Including this effect, a new corrected value is 
derived, which we display as a dashed line in Figure [4] 

However, for model A, which has a scaleheight h — 10 5 £0, 
a direct application of equations ( 17 1 or ([8]) yields a drag force 
Fa ~ 156 po Vq £q, which is more than four times the com- 
puted value. The main reason for this discrepancy is that the 
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FIG. 2. — Cuts of density and velocity along the z = plane after an integration time of t = 55§o/^0 for all the models (as indicated in the top labels). The 
gray-scale (color-scale in the online version) shows density, in units of the midplane density po, in logarithmic scale (as indicated by the bars on the right). The 
arrows represent the velocity field. The panels only show an inner region of the entire computational domain where the wake is formed. The perturber is located 
at (0, 0) and the axes are given in units of go- 



model considers all the mass in the disk, while the numerical 
simulation is restricted to the drag by streamlines with impact 
parameter less than L y = 7.5£o- In order to apply our ana- 
lytical approach to our simulation (or to a truncated Gaussian 
disk), the integral in the right-hand side of Equation |9]l should 
be performed between the lower limit I = 1 and the maximum 
impact parameter, in our case Z max = L y / (2£o) = 3.25. For 
model A it results in a drag force ~ 27 po Vq£q, which is 
the value plotted in Fig.^a) as a dashed line. For the remain- 
ing models B, C and D, this correction is marginal and the 
difference between the dashed, and the solid or dotted lines is 



due to the correction in the accretion rate. 

In general we see that the agreement between the analytical 
model and the numerical values is good. The results for the 
higher Mach number model (D) agree slightly better, which is 
expected since the ballistic model is applicable to hypersonic 
flows. 

4. DISCUSSION: APPLICABILITY OF THE MODEL TO 
A PLANET INSIDE A KEPLERIAN DISK 

The models presented in this paper (analytical and numeri- 
cal) consider a perturber traveling on a straight-line trajectory 



Gravitational drag from a disk 



Model A : h=10 5 f (thick disk), Moch 5 



30 



20 



10 



-10 



^KKS>(!(WS>iK!iK>l>!*»KI<><><>l>ll<K! 



2 4 6 

Model C : h = 0.5f (thin disk), Moch 5 



5 H 



> 3 

o 

~ 2 

X) 



-1 



<x><> 



w ><><><><><,<><>< >( | 



2 4 
f tfo] 



Model B : h = 2f (thick disk), Moch 5 
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Model D : h = 0.5f (thin disk), Moch 10 
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FIG. 3. — Mass flux M (in units of po vq £q) as a function of the cylindrical radius £ of the control volume, in units of £o- The symbols represent the average 
value (at each £) of the last five outputs in the simulations, corresponding to evolutionary times from 51 t o 55 £0/^0- The error bars show the maxi mum departure 
from the mean value. The solid horizontal lines are the prediction from the analytical model: equation |10l for panels (a) and (b), and equation [20) for panels 
(c) and (d). The horizontal dotted line (overlaid with the solid line for A) shows the full solution from Eq! Jvl. The dashed horizontal lines show the full solution, 
but including a correction due to the limited extent of the computational box (see text for details). The vertical dotted lines depict the size of the perturber. 



at constant velocity within a medium with a Gaussian density 
stratification. The medium has been regarded as a disk as one 
of the motivations was to obtain the drag force on an object 
orbiting in such a system (e.g. planets embedded in proto- 
planetary disks, or black holes in galactic nuclei). In a real 
situation, the perturbers move in complex orbits, the density 
of the disk and the relative velocity of perturber with respect 
to the ambient gas change with position, and the flow is not 
plane-parallel but presents differential rotation (shear). One 
can use the 'local' approximation, that is estimating the drag 
force at the present location of the perturber ignoring the cur- 
vature of the orbital motion, density gradients and shear ef- 
fects, as long as the curvature radius is much larger than the 
scales involved in the model (i.e. £0 and h). 

The relevant scale for the drag force is h in the case of a 
gravitational perturber moving in a thick disk. Hence, since 
the relevant scale for curvature effects and shear is the size of 
the orbit, these effects can be ignored provided that the aspect 
ratio of the disk is small (see also Rein 2012). 

In the following we consider the case of a thin disk. Let us 
consider the example of a planet in orbit inside a Keplerian 
disk to check under which conditions the model is adequate. 
A parcel of gas in a Keplerian disk around a star of mass M* 
will have a velocity at a radius ro of 1 : 

/ GM, \ 1/2 
v k = ( ] . (26) 



locity v p and the disk is Vq = v 



r o J 

The relative velocity between the planet in orbit with a ve- 

1 This is only valid if there is no radial pressure gradient. 



and, therefore 

+ « 2 — 2v p ■ t>k = v 2 (a 2 + 1 — 2a cosp) 



where 



a = — 



(27) 



(28) 



and p is the local pitch angle (the angle between v p and the 
azimuthal direction). The planet is in corotation with the disk 
when p = and a = 1, 

Combining equations ([3j, ( p6| and (27 1, we can write the 
ratio between the orbital radius and the planet gravitational 
radius as 



£0 M, 



(a 2 + 1 — 2a cosp) 



(29) 



where M p is the mass of the planet. Since the large-scale 
effects such as shear and curvature of both the orbit and the 



disk will be important when ro < £o> we can use equation ( 29 1 
to obtain the space of parameters a and p for which £0 < 
and, hence, the local approximation is valid. This condition is 



-M* 



1 — 2a cos p) > 1 . 



(30) 



We see that when p is close to ir/2, the above condition is 
fulfilled for any value of a, provided that M* > M p , The 
most stringent condition for a occurs when p w 0, i.e. when 
the planet is close to pericenter and apocenter. Taking p«0, 
we find that the condition ( |30"| i is not met for a's in the interval 
between a_ to a+, where a± w M p /M*. If we take for 
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FIG. 4. — Fjj (in units of po Vq §q) as a function of the cylindrical radius £ of the control volume, in units of go- For large values of £, it corresponds to the 
total drag onto the perturber. The symbols represent the average value (at each £) of the last five outputs in the simulations, corresponding to evolutionary times 
from 51 to 55 ^o/^0. the error bars show the maximum departure from the value. The solid horizontal lines are the predic tion from the analytical mode, equation 
I for A and B (For A the value is ~ 156, which lies outside the range plotted, see text for more details) and equation |22j for B and C. The horizontal dotted 
i (also ~ 156 for A) shows the full solution from Eq. (9). The dashed lines include the additional correction for the finite extent of the computational box, see 
text for details. The vertical dotted lines show the size of tne perturber. 



instance M*/M p = 1000 (similar to the Sun/Jupiter system), 
it is simple to see that for orbits with eccentricities larger than 
0.2, it holds that a > a + at pericenter, whereas a < a_ 
at apocenter. Hence, we conclude that large-scale effects are 
small at any radius as long as e > 0.2. 

A similar analysis can be done using the Hill radius of the 
planet in orbit : 



r H « a (1 



V 3M* 



1/3 



(31) 



where a is the semi-major axis of the planetary orbit, and e 
its eccentricity. The Hill radius roughly defines the region 
around the planet in which its gravity dominates. Outside a 
sphere with radius rjj the gravitational pull from the central 
star is not negligible. Thus, in order to apply our model to the 
wake of a planet in a disk, the Hill radius has to be larger than 
a few times £o (the size of the wake). Replacing equation (29 1 
into ( (31} we can write 

2/3 

■ (32) 



Co 



a (1 - e) (a 2 + 1 - 2a cosp) / M* 



3V3 r 



The gravitational pull from the star cannot be neglected if 
r H Co- Again, the stringest constraint for a occurs when 
p ~ 0. In that case, the interval of a for which rn < Co, is 
a'_ < a < a'i, where 



c4 = 1 ± 3 1 / 6 



Lo(l-e) 



1/2 



AL 



1/3 



(33) 



, that is ro/(a(l — 
e)/(l — e) at apoas- 



The factor in the square brackets in Eq. (|33 
e)), varies between 1 at periastron to (1 
tron. Let us consider the example of M + /M p = 1000 in a 
highly eccentric orbit with e = 0.8. At the apoastron, we 
need a < a'_ = 0.64 to have £o < rn- This condition for a 

is fulfilled because v p /v\ = -^/(l — e)/(l + e) at apoastron. 
In general, £o < rn at any radius of the orbit provided that 
e > 0.2. In conclusion, we have shown that, for eccentric 
orbits with e > 0.2 or retrograde disks, the curvature effects, 
as well as the pull from the central star, on the wake can be 
neglected. 

5. SUMMARY 

We present an analytical model for predicting the accretion 
rate and gravitational drag on a point mass that travels hy- 
personically along the midplane of a stratified medium with 
a Gaussian vertical density profile. The model considers that 
the trajectories of fluid parcels are ballistic, and is a direct ex- 
tension of the uniform density case presented in Canto et al. 
(20111. The analytic model was then compared with a set 



of three-dimensional, isothermal hydrodynamic simulations. 
In contrast to the geometry used in Canto et al. (2011 1, in 
the present work cylindrical control volumes are used, which 
makes the treatment easier and more natural given the sym- 
metry of the problem. 
The results can be summarized as follows : 

(i) Fully analytic expressions for the mass accretion rate 
M acc and the gravitational drag force are obtained 
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assuming a free streaming ballistic flow. 

(ii) Simpler expressions for M acc and are obtained for 
two asymptotic limits; for a thin disk, in which the disk 
scaleheight is small compared to the gravitational radius, 
and for a thick disk, in which the disk scaleheight is large 
compared to the gravitational radius. 

(iii) In a thick Gaussian atmosphere with vertical scaleheight 
h, and surface density £ the drag force, including the 
contribution of the nonlinear part of the wake, is given 
by 



F, 



d, thick 



4 v / 7r(GM) 2 E 



hvl 



In 



y/eGM/(2v%) and r n 



where r m \ 

while the mass accretion rate is 

A^{GM) 2 Y 



acc, thick 



hvl 



(34) 



e 2 h/{2^), 



(35) 



(iv) We find that the mass accretion onto a point mass M 
moving in a thin disk is proportional to M and to the 
surface density Y of the disk, 



M, 



acc, thin 



AG MY, 



(36) 



(v) The gravitational drag force in the thin disk regime is 
independent of the flow velocity (provided that the flow 



is hypersonic), and is also proportional to MY. 

-Fd.thin ~ 2ttG MY . (37) 

At the same time the gravitational deceleration of the hy- 
personic perturber in a thin disk (Fd/M) is independent 
of its parameters (i.e. mass or velocity) and depends only 
on the surface mass density of the disk. 

(vi) The mass accretion in the simulations was highly vari- 
able. The analytic prediction is slightly higher than 
the average values obtained from the simulations, but 
well within their dispersion. These differences might be 
due to the relatively low resolution of our present three- 
dimensional simulations. 

(vii) The drag force computed from the simulations showed 
a good agreement with the analytic model, with a lower 
dispersion compared with the accretion rates. The result 
(from the analytic model) that, in a thin disk, the drag 
force is independent of the Mach number of the flow is 
confirmed by the numerical simulations. 
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